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A mathematical model which describes the distribution of the temperature field and stresses originated upon 
cooling of a glass tube is proposed. The model is developed proceeding from the solution of a conjugated heat 
transfer problem The model is presented in the form of a system of finite-difference equations and makes it 
possible to analyze the dynamics of the annealing process in various cooling regimes using computer simulation. 


The process of annealing that determines the character of 
distribution of temporal and residual stresses in a glass arti¬ 
cle plays an important role in glass tube production. Unsatis¬ 
factory distribution of residual stresses reduces the service 
characteristics of the glass tubes and results in their rejection. 
Exceeding the rated thermoelastic strains can lead to sponta¬ 
neous fracture of the article in the stage of annealing. A 
mathematical model should be formulated to substantiate the 
choice of the annealing regime. The models developed ear¬ 
lier provide for determination of the temperature field and 
stresses in the wall of the glass tube (or any other article of 
the same shape) upon static cooling [1-3] they did not al¬ 
low describing the dynamic aspects of the heat transfer pro¬ 
cesses and initiation of stresses. 

Thus, the goal is to develop a mathematical model which 
describes the heat transfer processes, initiation and distribu¬ 
tion of stresses in a glass tube wall upon annealing. The 
model is also adaptable to optimization of annealing of glass 
tubes and other cylindrical glasswork, e.g., bottles. 

Consider a glass tube of radii R x and R 2 (Fig. 1). Let the 
temperature gradient be across the tube wall alone. Thus the 
temperature field obeys the equation 

ar,(r,0 X, fd 2 7j(r,/) j dT x (r,t)) 

- -=—I--— +--- , (1) 
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where 7j (r, f) is the temperature of the glass layer; r is the 
layer radius ( R 2 <r< /?,); t is the current time; X, is the heat 
conduction; c, and p, are the heat capacity and density of 
the glass, respectively. 

The initial conditions are: 

T\ ( r * *o ) - V 

where t Q and T 0 are the initial time and initial temperature 
of the glass, respectively. 
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The boundary conditions outside the tube obey the equa¬ 
tion 

ST (r,f) / x 

‘■(r.t'.o-r.to). (2) 

where r = R x is the external diameter of the tube; is the 
heat transfer coefficient on the exterior surface of the tube; 
and T (f) is the ambient temperature. 

The boundary conditions inside the tube are given by the 
following system of equations: 

, 57> . ST 2 (r,t) 

■'~dr~ 2 ~dr'’ ( 3 ) 

r, (r,/) = r 2 (r,0 

where r = R 2 is the internal diameter of the tube; X 2 and 
T 2 (r, f) are the heat conduction and air temperature inside the 
tube, respectively. 

The analysis of boundary conditions (3) leads to the ne¬ 
cessity of calculating the temperature field of the air inside 
the tube which obeys the equation: 

ST (r,f ) X (d 7 T 2 (r,i ) , 8T 2 (r.f )) 

-- - --— -r- + -- - , (4) 

dt c 2 p 2 ^ dr 2 r dr 

where 0 <r<.R 2 ; c 2 and p 2 are the heat capacity and den¬ 
sity of the air. 

The initial conditions are given by the equation: 

Tj / o) _ T o ■ 

The presence of complicated conjugating boundary con¬ 
ditions (3), uncommon to heat transfer theory, is attributed to 
the fact that air-glass heat transfer in the closed space of the 
tube proceeds mainly through heat transmission rather than 
via convection because the temperature is constant through¬ 
out the inner surface of the tube (with no regard for heat 
transmission at the ends of the tube). 
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Let us divide the inner air space and the tube wall into N within the range of (0, R 2 ) and (7? ? , /?,) (see Fig. 1.). 
layers with spacing of Using network method [4] with due regard for boundary 

conditions (2), (3) we obtain from Eqs. (1), (4) a system of fi- 
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R- R< nite-difference equations for calculation of the temperature 
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field of the glass tube: 
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. j = 


T 2 y,k+\]=T { [i,k+\], j=l,i = N, 


where i and j are the sequence numbers of the layer (/', 
j = 1 ... N) (relevant to glass and air, respectively); k is the ceeding from the convergence of the solution of system (5): 


The calculation step (in time) should be determined pro¬ 


calculation step; T x [i. A:] is the temperature of the i-th glass 
layer at step t, k, (7", [i, k ]) is the effective heat conduction 
of the glass 

A.. (7\ [i, k]) = 1.08 +4.19 x 10 -4 


At a 


C 2 P2 ^2 
3X, 


+ 1.82 x 10 "' 6 (R x - R 2 )° 73 (r, [i, k ] - 273)“ 
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An algorithm for calculating the temporal and residual 
stresses in the glass tube wall was developed on the basis of 


calculated with allowance for the radiation component of the data obtained in [5 - 7], The same problem was solved 

heat transfer; T 2 [/, k] is the temperature of the y'-th air layer there for a glass plate using the relaxation theory of glass 

at step k; r, [j] = (R { - Ar, i) is the radius of the i-th glass melting. The algorithm developed is represented by a system 

layer, and r 2 [/] = (R 2 - A r 2 f) is the radius of they'-th layer of air. of finite-difference equations: 
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Ac[i, * + 1] = (Jr d [i, * + 1] - Ar [* + 1 ]) -r rf [/, k ] ); 
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Fig. 1 . Sectioning of a glass tube. 


where t/ [i, A] is the fictitious temperature of the i-th glass 
layer at step k, x , [z, k ], T , f [z, k ] ^ is the time of ficti¬ 
tious temperature relaxation / , s. 

vlTfrk ],7/[i,A ]) 

]J/M '• 


K, 


r{T,[i,k ],T/ [i,k ]) is the viscosity of the glass deter¬ 
mined from the set of equations: 

n(r,[iU ],t/ [i.k ]) = 


= A + 


B , 


T/[i,k ]-C 




T ,[i ,k ) T, f [i,k ] 




where A, fl,, B 2 ,C, K ] are constant for each glass compo¬ 
sition; b y is the index of fictitious temperature relaxation 
and depends on the glass composition; r d [z, k] is the 
dilatometric radius of the z-th glass layer at step kr, r Q [/] is 
the radius of the z-th glass layer at the beginning of the calcu¬ 
lation; a ( and are the temperature expansion coefficients of 
the glass in the plastic and glassy state, respectively; r,[z, A] is 
the free radius of the z-th glass layer at step k, Ary[i, A] is the 
change in the free radius of the z'-th glass layer at step k, 
A r [A:] is the actual change in tube radius at step A; a n [z, A] is 
the sum of stresses relaxed in the z-th glass layer at step k, 
Aa[z, A] are the stresses originated in the z-th glass layer at 
step A; E is Young’s modulus; p is Poisson’s coefficient; 
o[i. A] is the current value of the stress in the z-th glass layer 
at step A; x 2 ^r,[z',A ],t/ [z.A ]) is the stress relaxation 


time determined by 
t 2 (7’,[«',A ],j/ 



where K, is the constant which depends on the glass com¬ 
position; and b 2 is the index of stress relaxation which also 
depends on the glass composition. 

Figure 2 illustrates the process of glass tube annealing cal¬ 
culated with the following parameters of relaxation models 
(5), (6): N= 10, R x = 0.26 m, R 2 = 0.25 m, c t = 1200 J/(kg ■ K), 
p, = 2500 kg/m 3 , c 2 = 1000 J/(kg • K), p 2 = 0.5 kg/m 3 , Aj = 
0.06 W/(m • K), E = 1 • 10 10 Pa, p = 022, a, - 3.3 • 10 " 5 K ~ 1 , 
02 = 1.1 • 10- 5 K-‘, h, = 0.65, b 2 = 0.55, AT, = 10 98 , 


t,°C 

550 


450 


350 


250 


150 


50 


0 500 1000 1500 2000 2500 3000 3500 4000 4500 t, sec 

a, MPa 

28 
21 
14 
7 
0 
-7 
-14 
-21 
-28 

0 500 1000 1500 2000 2500 3000 3500 4000 4500 x, sec 

Fig. 2. Simulation of the glass tube annealing process: J) air tem¬ 
perature in the center of the tube (along the cylinder axis); 2 ) tem¬ 
perature of the inner surface of the glass tube; 5) air temperature out¬ 
side the tube (ambient temperature); 4) stresses in the external tube 
surface; 5) stresses in the internal tube surface. 

A =-2.78, B x = 5574K, fl 2 = 15,500K, C=474K, 7o = 873 K, 
A, = 50 W/(m 2 - K). 

The solution of the problem considered makes it possible 
to develop a mathematical model for describing the process 
of stress initiation and distribution of temperatures and 
stresses upon cooling of a glass tube. The model proposed 
makes it possible to analyze the dynamics of annealing in 
various cooling regimes using computer simulation, thus 
providing for an evaluation of the efficiency of the annealing 
regime. The model can be used to analyze the dynamics of 
annealing in different cooling regimes based on computer 
simulation, making it possible to assess the efficiency of an 
existing annealing regime or comparing and then selecting 
the most efficient regime (process control). The model devel¬ 
oped can also give grounds for optimal process control, i.e., 
for optimization of cooling on the basis of selected criterion 
of quality within prescribed limits. 
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